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ABSTRACT 

In this paper the Yennie-Frautschi-Suura approach is used to simulate real 
and virtual QED corrections in particle decays. It makes use of the uni- 
versal structure of soft photon corrections to resum the leading logarithmic 
QED corrections to all orders, and it allows a systematic correction of this 
approximate result to exact fixed order results from perturbation theory. 
The approach has been implemented as a Monte Carlo algorithm, which a 
posteriori modifies decay matrix elements through the emission of varying 
numbers of photons. The corresponding computer code is incorporated into 
the Sherpa event generator framework. 
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1 Introduction 



The next round of collider-based experiments in particle physics will start with the LHC finally 
becoming fully operational, and producing particle collision data at unprecedented rate and 
energy. In the preparation for this huge enterprise, a new generation of Monte Carlo simula- 
tion tools, hke Pythias p], HerwiG++ [2] and Sherpa [3] has been constructed, to meet the 
increasingly complicated experimental situation and the demand for an improved description 
of data at a higher level of accuracy. In addition, it was anticipated that by moving to the new 
programming paradigm of object orientation and modularization the more mundane software 
management task of code validation and maintenance could be addressed in a more transpar- 
ent and alleviated way. This leads to the typical strategy of event generators, to dissect the 
simulation of full events into different physics aspects, being better reflected in the modular 
structure of the emerging new codes. 

In this paper the construction of a new physics module for the Sherpa framework is dis- 
cussed, which deals with the simulation of QED radiation in particle decays. Up to now, this 
has typically been left to the PHOTOS [IfS] programme. However, there have been two reasons 
for replacing PHOTOS: First of all, PHOTOS builds on a parton-shower like coUinear approxi- 
mation for the simulation of photon emissions, which intrinsically has some shortcomings when 
the mass of the decaying particle is not much larger than the masses of its decay products. This 
has already been noted in pi7] and triggered the development of the module SOPHTY [6J in the 
framework of the HERWIG++ event generator. It also has become a wide-spread belief among 
the authors of the main event generators that the maintenance of interface structures to addi- 
tional codes such as PHOTOS, supplementing QED radiation to the simulation, overwhelms the 
burden of constructing and maintaining corresponding modules directly in the event generators. 

Similar to the case of SOPHTY in HERWIG++, the construction of the new module, PHOTONS++, 
in Sherpa bases on the approach of Yennie, Frautschi and Suura (YFS) p] for the calculation 
of higher order QED corrections to arbitrary processes. This approach resides on the idea of re- 
summing the leading soft logarithms to all orders, rather than focusing on the leading coUinear 
terms. These soft logarithms are largely independent of the inner process characteristics and 
can be calculated from the external particles and their four-momenta only. The big advan- 
tage of the YFS formalism is that in addition it allows for a systematic improvement of this 
eikonal approximation, order- by-order in the QED coupling constant. This explains why a good 
fraction of the most precise tools for the simulation of QED radiation root in this algorithm 
[9,10fTT][T^. 

In contrast to the SOPHTY implementation the aim of this implementation is to address also 
decays with more than two flnal state particles. This leads to different strategies of enforc- 
ing four-momentum conservation after the soft photons are reconstructed. In addition, some 
correction terms to restore precise results for the flrst order in the electromagnetic coupling 
constant are employed, improving the formal accuracy of the results of PHOTONS++. This also 
has not been included in SOPHTY. 

The outline of this paper is as follows: After briefly reviewing the YFS formalism in Sec. [2] 
in the framework of particle decays, the Monte Carlo algorithm adopted here is detailed in 
Sec. El Then, some higher order corrections are discussed in Sec. IH Finally, in Sec. EJthe new 
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code is validated through a detailed comparison with HORACE and WiNDEC [13] for the case of 
leptonic Z and W decays before some results relating to other particle decays are presented. 



2 YFS-Exponentiation 

In this section, the YFS approach [8] for an approximative description of real and virtual 
QED corrections to arbitrary scattering or decay processes will shortly be reviewed, in the 
framework of particle decays. The virtue of this formalism is that it can systematically be 
improved, order by order in the electromagnetic coupling constant a. The YFS approach bases 
on the observation that the soft limits for matrix elements with real and/or virtual photons 
exhibit a universal behaviour, and on the fact that the corresponding soft divergences can be 
factorised into universal factors multiplying leading order matrix elements. 

When defining the final state as a configuration of primary decay products with momenta 
and any number of additional soft photons with momenta k the fully inclusive decay rate 
reads 



1 °° 1 /" 



oo 

E -^"r^"" 

nv=0 



where pi is the four-momentum of the decaying particle. Here and in the following ny and riR 
are the numbers of additional virtual and real photons, respectively, that show up in the higher- 
order matrix element but not in the uncorrected zeroth order matrix element (thus labelled by 
A^o). The starting point of the YFS algorithm is to approximate these dressed matrix elements 
through the zeroth order one times eikonal factors, which depend on the external particles only. 
This effectively catches the leading logarithmic QED corrections to the process. The correct 
result can then be restored order by order in perturbation theory by supplementing the non- 
leading, process-dependent pieces. 

In the case of one virtual photon this can be formalised as 

Ml = aBM^o + , (2) 

where Mq is the infrared-subtracted matrix element including one virtual photon (with Mq 
finite when k —>■ due to the subtraction). All soft divergences due to this virtual photon are 
contained in the process- independent, universal factor B, see Appendix El for a more thorough 
discussion. Here, and in the following, the sub- and superscripts denote the number of real 
photons and the order of a, respectively, both for the infrared-subtracted and for the original 
matrix elements. 

Similar to the one-photon case, YFS showed that the subsequent insertion of further virtual 
photons in all possible ways leads to 



= Mo" 
Ml = aBM^ + Ml 

Ml = ^^^M',+aBMl + Mi 



(3) 
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and so on. Therefore, for a fixed order in a 

nv 



ny—r 



JaB) 



r=0 



(4) 



and, summing over all numbers of virtual photons ny, 



nv=0 nv=0 



(5) 



Since photons do not carry any charge and because virtual photons inserted in closed charged 
loops do not produce any additional infrared singularitjfl, this can be generalised to any number 
of real photons, such that 



nv=0 



exp{2aB) 



nv=0 



riR 



(6) 



Hence, is free of soft singularities due to virtual photons but it still may contain those 

due to real photons. 

YFS showed in [S] that the factorisation for real photon emission proceeds on the level of 
the squared matrix elements rather than on the amplitude level. For a single photon emission 
therefore this yields 
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2(27r)= 



oo 


2 


oo 




= S{k) 




nv=0 




nv=0 



(7) 



nv=0 



Here, S{k) is an eikonal factor containing the soft divergence related to the real photon emission, 
see Appendix |X1 Denoting with P^^'^"'^ the complete IR-finite (subtracted) squared matrix 
element for the basic process plus the emission of ur photons including ny virtual photons and 
using the abbreviation 



nv=0 



n-v+riR 
riR 



(8) 



the squared matrix element for real emissions, summed over all possible virtual photon 



^ A similar program cannot directly be translated to QCD, where the emitted gluons act as parts of antennae 
emitting further gluons, thus modifying the pattern of possible infrared poles and thus leading logarithms in 
each emission. 
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corrections, can be written as 



2(27r)= 



i=l 

nR 

+E 



nv=0 

nR 



+E 



i=l 
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n 



nR 



PnR-l{kl, ■ ■ ■ , ki-i, ki^i, . . . , knji) S{ki) 



<>3 = 1 



P2{ki, kj 



nR 

n 



1=1 



S(h)S(kj) 

~l~ PnR{kl, ■ ■ ■ 1 knj^ 



~S{ki) 



+ ... 



(9) 



Demanding agreement with the exact result up to 0{a), this expression thus contains only 
terms with (3^, (31 and I3\. Then 



2(27r)s 



nR 



nv+2nR 
nR 



nv=0 
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(10) 



Inserting this into the expression for the decay rate and expressing the (5-functions ensuring 
four-momentum conservation as exponentials yields, 



2M-r = 



nR I 



k 



J d^y J d$p^ < exp [2aB] exp 



k 



X 



K 
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As before, all singularities due to virtual photons are contained in B, while all singularities 
due to real emissions are incorporated in the integral over S{k). To restore the momentum 
conserving (5-function the divergences have to be split off this integral. This can be done by 
simply subtracting the terms that are divergent for A; — > 0. To this end, a small "soft" region 
is defined together with an infrared-safe function Z}(f])0, such that 



^ s{k)e-'y' 



where 



and 



~S{k) 



2aB{n) + D{n) 



1 - 0(A;, Q) ] + e~'y''e{k, Q) + ( e-'^'' - 1 ) ( 1 - e(A;, Q) 



(12) 



-iyk 



- 1 



^S{k) 



^S(k) e-'y''Q(k,Q) 
k 



e{k,n) +e-'y''e{k,n) 



^s{k){i-e{k,n)) 



^S{k) 



(13) 



(14) 



Reinserting this into the cross section, executing the ^/-integration and reexpanding the exponen- 
tiated integral yields 



2Mr 



y— 

riR 



d$,,d$U2vr)V (Y,P^ - J2pf -J2k)e 



^2a{B+B{Q.)) 



riR 



I S(^ki 



(15) 



The whole factorisation is independent of possible spin correlations in the "undressed" matrix 
element. Thus, the same result is obtained if the spin-summed and averaged matrix element 
squared \A4f is replaced by Pai3.M°'.M^* where Pai3 is a spin density matrix. 
The infrared subtracted squared matrix elements read, up to 0{a), 



Pi 



1 



2(27r)^ 



1 * 

2 



(16) 



^ Obviously 8(fc,ri) divides the phase space into two regions. While comprises the region containing the 
infrared divergence, (1 — $7) is completely free of those divergences. Hence, Q{k,fl) — 1 if k ^ ft and zero 
otherwise. Thus, D{fl) is IR save and B{n) contains the divergence. 



7 



or 



/3° = p.,M°"M°^* 



k = P.p(^l^yrMf -~S{k)M^,-Mry (17) 

3 The Algorithm 
3.1 The master formula 

The basic, undressed matrix element (no additional photons) reads 

2M ■ To = j d<l>, {2'kY8\pc +Pn-Qc-Qn)\M\^ (18) 

where the differential phase-space element for the outgoing momenta q G {Qc, Qn} is given 
by 

Here, and in the following, the initial and final state momenta have been classified to whether the 
respective particles are charged or neutral: the sums of all initial state momenta are labelled 
by pc,N for charged and neutral particles, respectively, while Qc,n denotes the sums of all 
charged or neutral final state momenta. After QED corrections, the Qc and Qtv will become 
Pc and P/v, respectively. K is the sum of all additional real, resolved Bremsstrahlungs-photons 
generated in the process, whereas photons already present in the core process are included in 
P/v and Qtv, respectively (an example for this seemingly unlikely case would be the rare decay 
B+ K*+-f). 

In the previous section the factorisation of infrared divergent terms and the construction 
of infrared-finite expressions for cross sections with all possible numbers of resolved photons 
has been discussed. In these expressions the universal, process-independent parts of the QED 
corrections have been separated and exponentiated, the residual process dependence and the 
effect of particle spins etc. has been absorbed in infrared-finite, subtracted terms P, cf. Eq. (fT^ . 
With small changes in the notation this form of the cross section thus reads 

1 r 

2M-T = / d<l> e^(^) JJ5(A;i)e(A;„^])/3°C. (20) 

rij i=l 

Here, the phase space has been separated into a phase space element for the particles of the 
"core" process and one for the additional resolved real photons, 

d<l> = d% d<Dfc (271)^5 (pc +Pn-Pc-Pn-K). (21) 
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with 

d*. = n^- (23) 

i=l 

Note that the factor 2(2^' missing in the photon phase space element, has already been incor- 
porated in the eikonal factor S{k), in accordance with the choice made in [8]. In the equation 
above, Eq. ( l20i) . the undressed matrix element /3g has been factored out and the remainder of 
the perturbative expansion in a has been combined in the factor C, 

Furthermore, the YFS-Form-Factor has been introduced 

Y{Q) = y^A^) = 5Z 2a (5,, + 5,,(f])) (25) 

i<j i<j 

where the sum i < j runs over all pairs of charged particles, taking into account each pair only 
once. The infrared factors Bij and Bij are defined as 

- j d k- (^^, _ . ^^^^ + ^2 + 2(A: . p,)9j ^^^^ 

B,,{n) = -l^z,z,e,e, J d^kS{e){i-e{k,n))(^-^-^^ . (27) 

They are the generalisation of the quantities defined in the last section, cf. Eqs. ([2]) and (IT^ . 
Both contain the virtual and real infrared divergences, respectively. These divergences cancel 
according to the Kinoshita-Lee-Nauenberg theorem [T^fT^ . Thus, each Yij{Q) is guaranteed 
to be finite, which is explicitely shown in Appendix |X1 In the terms above, Zi and Zj are the 
charges of the particles i and j in terms of the positron charge e, and the signature factors 
6 = ±1 for particles in the final or initial state, respectively. The symbol G, already defined at 
the end of section [21 refers to a phase space constraint with Q denoting the soft, unresolvable 
region of photon radiation. Hence, Q{k,Q) = 1 if k ^ Q and zero otherwise. If this division is 
done by defining an energy cut-off, the definition of Q is not Lorentz-invariant and the frame in 
which this cut-off forms a flat hypersurface also needs to be specified. The advantage of splitting 
the photon phase space in that manner lies in the alleviation of integrating S{k) over k. If the 
cut-off is defined in the frame the photon generation and momentum reconstruction will be 
done irj^l then the integration over the photon energy separates from the angular integration 
(see Appendix [nj , leading to yet another simplification of the calculation. 



In the algorithm presented here, this will be the rest frame of the multipole, i.e. the combined rest frame 
of all charged particles pc + Pc ■ 
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s{k) = > ; s,,{k) = -^y: z,z,e,e, ( - A: ) • (2^ 



The eikonal factor S{k) has aheady been introduced in the last section. It is defined as 

i<j i<j ^-^^ 

However, despite all terms being finite in Eq. (1201) . it cannot be used straight away for 
Monte Carlo generation. This is because it is written in terms of the already corrected final 
state momenta pi and not the original undressed momenta g^. The problem here is that the 
undressed momenta are defined in an n-body phase space whereas the dressed momenta are part 
of an {n + ?T,^)-body phase space. This neccessitates a mapping procedure of the n-body onto 
the (n + n^)-body phase space. In principle, details of this mapping procedure are irrelevant 
as long as it respects the soft limit of photon radiation not altering the original kinematics, i.e. 
in this limit the momenta of the orinial particles in the (n + n^)-body phase space have to fall 
exactly onto those of the n-body phase space. 



3.2 Phase space transformation 

To solve this, consider the rest frame of all charged particles involved in the basic matrix element 

Pm=Pc + Pc- (29) 

These particles form the multipole responsible for the Bremsstrahlung of the additional photons. 
In the rest frame of this multipole, a simple form of the mapping consists of a mere rescaling 
of the three-momenta of all final state particles by a common factor u such that the additional 
photons are accomodated. Clearly, the initial state momenta cannot be altered, because they 
have already been fixed when the basic matrix element was calculated. So, the task is to 
rewrite Eq. fl2Tl) . explicitely in the rest frame of the multipole in question. The neccessary 
transformations are detailed in the appendix, cf. App. [BTTI here it suffices to give the result. It 
reads 

d* = d<l.,d<l.i,(2,r)*^(po+PA'--Po--PA'- A') 



n 



(2vr)32p0 



n 



(27r)^5 {Pc+Pn-Pc-Pn-K) 



d^p d$i 



m 



M,p 



{2'nf5\P^,) {271)5 {P\ 





M 



pO 



(30) 



In a similar fashion, the phase space related to the zeroth order uncorrected cross section can 
be transformed to 

d$o 



(2vr) 



M$g 5^ {pc +Pn-Qc- Qn) 



m 



M,q 



■ d% {27rf 6'{Qm) {27r)6 (Q°, - Q°; - p'c) 



In both cases, nijv/,p {^M,q) is the invariant mass of the corrected (uncorrected) multipole and 
the vector components P^ and P^ {Q^ and Q^) are taken in the Pm (Qm) rest frame. The 
Jacobian emerging in both cases will ultimately find its way into a correction weight in the 
Monte Carlo realisation of the method. 
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3.3 Mapping of momenta 



As mentioned before, the mapping procedure still has to be defined in detail to reconstruct the 
particles' momenta. The basic ideas of the mapping procedure suggested here are as follows: 
When representing all four-vectors in the rest frame of the multipole 

• treat all final state momenta equally 

• scale their three-momenta by a common factor u 

• distribute the photon momenta 

• assign the energy- component of every vector such that momentum conservation and all 
on-shell conditions are fuUfiUed 

This will ultimately necessitate a change of the initial state momenta as well. However, since 
they are already fixed for the calculation of the basic matrix element this change will reduce to 
employing another frame during the reconstruction procedure. 

However, closer examination reveals that the mapping paradigm above in fact enforces a 
different treatment for purely neutral and partially or fully charged initial state configurations. 
The reason is that the momenta of the newly generated Bremsstrahlungs photons need to 
be balanced. Furthermore, the phase space integral still has to be rewritten in terms of the 
undressed, original final state momenta defining the original matrix element and cross section 
without QED radiation. This will be addressed in the next sections. Sec. I3.3.1ll33^ where 
the case of decays, i.e. single initial state particles, either neutral or charged, will be discussed 
separately. Formally, of course, both treatments will yield identical results, since only the soft 
limit of photon emission is defined from first principles and because both treatments respect 
this limit. 

3.3.1 Neutral initial states: final state multipoles 

The first case to be considered is the case of a neutral particle of mass M decaying into a 
final state with charged particles. The reconstruction paradigm above completely fixes the 
reconstruction procedure to a rescaling of all final state momenta, both charged and neutral, 
and balancing the summed photon momentum K by moving the frame of the multipole and, 
hence, of the initial stat^. Denoting, again, with the undressed and with pi the dressed final 
state momenta, and denoting their respective sums by Qc, Qn, Pc and Pn, as declared earlier, 
and using K as the summed momentum of all Bremsstrahlungs photons, the reconstruction 
prescription reads as follows: 

Note that it is not possible to distribute any fraction of the photon momentum equally to all final states 
with the constraint that the multipole remains in its rest frame. It therefore is mandatory to balance the photon 
momentum with the initial state. 
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• The momenta of the Qm rest frame 

Q% = {Q%.Qn) (32) 

will be mapped onto 



Pn 



p'n' = 


M2 + {uQn + uQn + K^upN + K^ 




— 




(33) 






(34) 


K" = 




(35) 



in the Pm rest frame. 

Pn and p'^y are the same physical vector but in different frames. The scaling parameter u 
now is determined by momentum conservation, i.e. 



C N 

where the subscripts C and N in the sums indicate a summation over charged and neutral 
particles, respectively. 

• The phase space element expressed in terms of the undressed final state momenta then 
reads 

3 



d$ = {27r)U%d^k S'{QM)s{Q'i,-Q''c-p'h) 

(37) 



M2 (po + po + K^) 

n r 

Qi 



T,C,N^ i-- 

3.3.2 Charged initial state: Mixed multipoles 

The other case of relevance in the framework of this publication is the decay of a charged particle 
of mass M, leading to multipoles containing both initial and final state particles emitting the 
photons. Again the paradigm above completely fixes the reconstruction procedure. Basically, 
the problem is to compensate the photon momentum after the final state momenta have been 
rescaled. This is achieved in the following way: 



pI 
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• The momenta of the Qm rest frame 



will be mapped onto 



Pc ^ Pc 



(38) 



+ {uQc — nc^y, —uQc + ncK, = upc + ncK, 
Pq, uQc - ncfij 

(^P^, uQm - riNK = -2uQc - unk^ 

i?) (39) 



in the Pm rest frame. Here, nc and un are the numbers of charged and neutral final state 
particles, respectively, and the abbreviation 



K 



2nc + un 



-K 



(40) 



has been introduced for a more compact notation. Again, pc and p'q are the same physical 
vector represented in different frames, thus specifying the relation between the Qm and 
the Pm rest frame. In the soft limit, i.e. for X — > 0, the scaling parameter — > 1 and 
both vectors are identical, as required. 

• In general, the scaling parameter is fixed through energy conservation as the solution of 



= ^J^pV(u^, 



C N 

(41) 

The phase space integral rewritten in terms of the Qi reads 

d$ = {2nr d$, d$. 6' (Qm) S {Ql, -Ql- p^c) (po +^p^ + ^0) 



X u 



3n-4 ^ ^C,N ^ 



■PC _ rMi. 11- 

7^ 2^C,N pO i=l 



Pc 



Pi 



(42) 



It is worth noting that this is identical to the case of a neutral particle in the initial state. 
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3.4 Event generation 

Having transformed the phase space integrals allows to write the full decay rate including real 
and virtual QED radiation as 

2^ ■ r = E ^TT / "^"^^ dM^rryS' (Qm) 6 {QIj - " Pc) ^'^^''^ Po C 



n 

i=l 



n r 



n 



pY i=l 



Pi 



(43) 



where p and p' now generally stand for the initial state particle. 

The zeroth order differential decay rate dFo, which will be used by default in all decays in 
Sherpa can easily be extracted and, employing Eq. ( pTl) . reads 



r = 5^— / drod<i>,e^(^)n p(^^)®(^-^) 

' 1=1 



m 



M,p 



Q'c + Q^ 



u 



3n-4 P° ^ 



C,N ^ 



n 



C,N ^ i=l 



.Pi 



c. 



(44) 



Up to here no approximations have been made at all. In order to generate the corresponding 
distribution with Monte Carlo techniques, however, this form is not particularly useful. To 
simplify Eq. (jUj) therefor, hit-or-miss and reweighting techniques are used, demanding upper 
bounds for the various pieces: 

• All higher orders are neglected, thus setting C to one. 

• The maximum of all Jacobians is given for K = 0, coinciding with the leading-order phase 
space. 

• The dependences on the dressed momenta in the eikonal factors are removed by approxi- 
mating these factors by those depending on the undressed variables from the generation. 



The resulting crude distribution reads 



n^=0 i=l 



(45) 



After executing all ^-integrations giving 



n^4(A;,)e(A;„fi)=n"- 



(46) 



i=l 



the YFS-Form-Factor is estimated by 



-n 



(47) 
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for suitable choices of f2 0. Reinserting this into the crude estimate, the leading-order cross 
section can be seperated from the QED radiation, and 



-e "n"-^ 



n-y=0 



(4J 



The result is the undressed zeroth order cross section times a Poisson distribution with the 
avarage photon multiplicity n. In this factorised state the photon distribution can be separated 
from the generation of the basic matrix element. Assuming the latter to be already generated it 
can a posteriori be corrected to the leading-logarithmic all-order QED correction by generating 
the photon distribution as follows: 

1. Generate the number of photons according to a Poissonian distribution with mean n. 

2. Generate each photon's momentum according to Sq{k). This implies 

• that its energy is distributed according to 

P{k') - ^ (49) 

• and that the azimutal and polar angles are distributed according to 

p{9,c^)^Y.(— ' (50) 

where is a null vector of unit length, 

e^ = -^F with el = 0. (51) 

It is possible that more than one hard photon is created such that the total energy of 
all photons exceeds the decaying system's energy. Obviously, this has to be avoided to 
guarantee energy conservation. A simple way of achieving this is a mere veto on such 
situations, accompanied with a repetition of photon generation, starting from step 1. 

3. Reconstruct the momenta. 

4. Calculate and apply all weights. This yields a total weight, namely 

W = l^dipole X Wyps X Wj,l X Wj,m X Wc , (52) 

In this publication (and in the code), this choice has been to hniit the photon energies by setting an infrared 
energy cut-off of O.lGeV, unless otherwise stated. 
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where the individual weights are given by 



W^dipoie = 11-- — -— (53) 

Wyfs = exp{Y{pc,Pc,^)+n) (54) 

p'O Z^CAf ^ i=l ^^-^ / 

l^c = C. (57) 

Here, PFdipoie corrects the emitting dipoles from the unmapped to the mapped momenta, 
Wyfs accounts for the exact YFS form factor, Wj^l essentially denotes the Jacobian due 
to the Lorentz-transformation, Wj^m is the weight due the momenta-mapping, and Wc 
incorporates higher-order corrections, where available. 

The maximum of the combined weight indeed is smaller than the maximal weight em- 
ployed for generating the distribution, W < W{K = 0). Hence application of the com- 
bined weight is just a realisation of a hit-or-miss method. The distribution obtained is 
now the exact distribution of fl20l) or (HH) . 



4 Higher Order Corrections 

In the last section, the procedure generating the QED corrections to cross sections, following 
Eq. f|T5l) . has been outlined. By construction, the algorithm yields exact all-orders results, if all 
matrix elements are known. This, however, is never the case. On the other hand, the dominant 
universal soft photon contributions, real and virtual, are included to all orders in the YFS 
form factor, Eq. ( !25|) . Thus, if the zeroth order undressed matrix element only is known, i.e. 
if C = 1, the photons will be solely generated according to a product of eikonal factors S{ki). 
Consequently, their distribution will be correct in the soft limit only. Away from this limit, 
exact matrix elements at a given order may be mandatory to yield satisfactory and sufficient 
accuracy. For most applications on decay matrix elements - the topic of this publication - it 
will be sufficient to implement the matrix element correction to the first order in a, i.e. for the 
emission of one additional real or virtual photon. It should be noted here that hard photon 
emission predominantly occurs in situations where potential emitters are characterised by a 
large energy-to-mass ratio and that in any case hard photon emissions tend to populate regions 
in phase space that are coUinear w.r.t. the emitters. In contrast, large angle radiation has the 
tendency to be predominantly soft. 
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4.1 Approximations for real emission matrix elements 

As already explained, the vast majority of hard photon emissions deserving an improved descrip- 
tion through corrections to the soft limit underlying the YFS approach occurs in the collinear 
region of emission phase space. In this region, the well-known collinear factorisation can be 
used to approximate f3l. This amounts to an inclusion of the leading collinear logarithms aris- 
ing in this limit, which are incorporated for instance in the Altarelli-Parisi evolution equation 
[TU] and corresponding splitting kernels. 

Since masses are to be taken fully into account the quasi-collinear limit defined in PTIIT^ 
replaces the more familiar collinear one. In this limit the matrix element factorises as 

, 2 r eX2g(-*)fe,fc)|A<0(p, + A;)|^ ifzGF.S. 
> \M?(pi,k) = < ^ (58) 

f I I e'Zfg(^-\p,,k)\Ml{x.p,)\' ifzGl.S.. 

Here the g(™'°"*)(p^j k) denote massive splitting functions. For instance, for the case of a fermion 
emitting a photon they are given by 



1 fp (,^_^ 

ip,-k) V ^ ip^-k) 

where x = '^o and z = ^^cq^ are the fractions of the fermion energies kept after the emission 
of the photon, and where Pff{y) is the well-known Altarelli-Parisi splitting function 

Pffiy) = Yzj- (61) 

The dipole splitting functions of [T7j have been generalised further in ^\ to incorporate 
also polarisation. Thus, in principle they could directly be used in the framework of the YFS 
formulation replacing the original eikonal factors. In the framework of this publication, however, 
they are employed as universal correction factors, reweigthing explicit photon emission such that 
the correct collinear limit is recovered. Since they interpolate smoothly between both limits 
they already include the soft limit. Therefore, in the correction weights, these soft terms have 
to be subtracted because they are already acounted for in the orginal YFS eikonals. In addition, 
since the dipole splitting kernels refer to an emitter and a spectator forming the dipole, for each 
dipole two such terms have to be applied, such that the squared matrix element with the dipole 
terms approximating the photon emission reads 

■^if = -e'J2[z^ZA^^g^M^PJ,k)\M'o\'] (62) 
' J2 [Z^ZA0J {giM^Pj, k) + g^,^p,,p., k)) \Ml\^] . (63) 
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Here, charge conservation in the form ^i^i — has been used. The second particle in each 
massive sphtting function g^j denotes the spectator of the emission process and accounts for 
the recoil, thus ensuring four-momentum conservation. It should also be noted that the sum in 
the equations above runs over charged particles only. 

In order to subtract the soft terms, it is useful to consider the soft and quasi-coUinear limits 
of the dipole splitting kernels gij{pi, pj, k): 



{Pi ■ k) V {Pi ■ k) + {pj ■ k) {pi ■ k) 



g.M^P.k) TT^ ( T, - TT^ ) (64) 



g,,{p.,p„k) ^--^ g(-*/-). (65) 

Because the soft limit is universal and spin-independent, it is a straightforward exercise to 
define soft-subtracted dipole splitting kernels 



gij {Pi ,Pj,k) = gij {puPj,k)- gff^^ {Pi,Pj,k) 



'^{Pi-Pj] 



m: 



2 



The soft-subtracted dipole splitting kernels g^j now have the correct (finite) soft limit while 
retaining the original quasi-coUinear limit of g^^ (Eq. fl65p ). Accordingly, the soft-subtracted 
matrix element can be approximated as 

= -£^Y1 ^^^A^j {gij{Pi,Pj, k) + gji{Pj,Pi, k)) (3^ . (67) 

i<j 

The exact form of the gij{pi,Pj, k) for different emitter-spectator configurations will be given 
in Appendix [Pi 



4.2 Exact real emission matrix elements 

In order to achieve an even higher precision, the implementation of exact higher-order full 
matrix elements becomes mandatory. It should be clear, however, that large differences with 
the approximated matrix elements above will occur only in non-singular regions of comparable 
hard, wide-angle emissions. Since the module presented in this publication, PHOTONS++, is 
embedded in the Sherpa framework it is easy to implement such infrared subtracted squared 
matrix elements, making use of tools and functions already provided within the framework. 
In particular, some basic building blocks for the calculation of helicity amplitudes already 
used in pUll^ can be recycled to construct the neccessary, infrared-subtracted one-photon real 
emission matrix elements, which are then evaluated at momentum configurations generated by 
the algorithm of Section [3l These building blocks are listed in App. [El Exact first-order matrix 
elements have so been implemented for a number of relevant matrix elements, see below. It is 
worthwile to stress that in principle also second-order precision could be achieved, if neccessary. 
In general, the infrared-subtracted squared matrix element can be written as 

Pi = ^T^.MfM!* -Sik)M'oM',* , (68) 
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and it is only the amplitudes M. that are process-specific and need to be listed for the different 
processes. It should be noted that within the Sherpa framework the real emission matrix ele- 
ments are straightforward to implement, in contrast, the incorporation of loop matrix elements 
is somewhat more involved: in those cases the integral has to be calculated analytically and 
the divergences must be cancelled before implementation as a function of the outer momenta. 



4.2.1 Two-body decays of type V — » FF 

The matrix elements for two body decays where one neutral vector particle decays into two 
charged fermions, V FF, reaqj 



Ml = ie e {p,X) u{pi,si) 



^ 7 rTV2 2^ ^^lPl + CrPr) 



[P2 + k) — 



v{qi,si) el' {k,K) . (69) 



Of course, momentum conservation must hold, and therefore p = gi + g2 in the former and 
p = Pi +P2 + k in the latter case. Hence, if of the generated event exceeds the number of real 
photons in the respective infrared subtracted squared matrix element, a projection of the higher 
dimensional phase space onto the lower dimensional one has to be performed. In practise, this 
amounts to redoing the reconstruction procedure using only a subset of all photons generated 
in that run. Furthermore, 

ClPl + CrPr = CL^^+CR^^. (70) 

Thus, the generic matrix element is adjustable to various decays of neutral vector bosons. A 
few key examples of the couplings cl and cr to the left and right-handed fermionic currents 
are listed in Table [H 

The real-emission matrix elements depend on the polarisations, and they are expressed in 
terms of the X, Y and Z functions listed in Appendix [El as 

Mo = ie X {qi,si]e^;q2,S2;cL,CR) (71) 



All particles involved are considered to be point-like, i.e. their vertices do not contain form factors. 
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Process 


cl 


cr 


Z 


2swcw ^ 






—ie 


—ie 



Table 1: Values of the coupling constants of different vector particles to the left- and right- 
handed leptonic currents. 



and 
Ml = 



le 



1 



+ 



2{pI - m?) 
1 



1 + 



m 



m 



2{pI - m?) 
1 

1 

2(pg - m2) 



'a J s 

X (pi, si, e:^'*, Pa, s, 1, 1) X [pa, s, e^,p2, S2, Cl, cr) 

s 

X (pi, Si, Pb, S, Cl, Cr) X {pb, S, S^*,P2, S2, 1, 1) 



^ ~ ] {P^^^'^^^^^Pb^^^CL,CR) X {pb,S,e"'*,P2,S2,l,l) 



(72) 



where 



Pa = Pi + k 
Pb = P2 + k, 



(73) 



and where pi and p2 are the momenta of the final state leptons. The bar over the fermion spin 
label Si signifies an anti-particle. 



4.3 Virtual emission correction /3q 

The only virtual corrections occuring to level 0{a) are 

For the above case of decays of the type V FF they read 



^0 



mi 



with 



A = 



1 in on-shell scheme 



I in MS scheme 



(74) 



(75) 



(76) 
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which agrees with [22|23] . Effects of potentially different left- and right-handed couplings cl and 
cr, cf. TabiH only enter in terms suppressed by ^ and are curently neglected in PHOTONS++. 
For the process W — > iu, cf. [21], the first order virtual correction reads 

^0 = f 

5 Results 

In this section some of the results of the PHOTONS++ module, as it is implemented within the 
Sherpa framework, are presented. The focus lies on the central distribution produced by the 
preceding calculations, the total energy of all photons radiated per event in the rest frame of 
the decaying particle. In addition, angular distributions for dipole and multipole configurations 
will be shown. 

5.1 Validation: leptonic heavy gauge boson decays 

The leptonic decays of W and Z bosons, W — > lui and Z II, will play the central role 
in validating the accuracy of the PHOTONS++ implementation of the YFS approach. Before 
studying in more detail these processes and comparing the results obtained with PHOTONS++ 
with those from other codes, namely AMEGIC++ [SH] and WiNDEC [21], it is worthwhile to 
discuss one of the key distributions, namely the total energy radiated off the decay. 

5.1.1 Radiated photon energy 

The result for this distribution, namely the total energy radiated off the decay in form of pho- 
tons, is presented in Figure [H For both processes, i.e. for both leptonic Z and W decays, 
different leptons with different masses have been considered. Clearly, radiation is most impor- 
tant in final states involving electrons, being the lightest fermions taken into account, while 
radiation off heavier fermions is increasingly suppressed. One of the most prominent features 
of every radiated energy spectrum is the kink at around half the boson mass, which is due to 
kinematics. It limits the energy involved in single photon emission off the final state fermion 
to its maximal energy, roughly half the boson mass. This kink gets washed out and moves to 
the left with increasing fermion mass. Events with total radiated energy surpassing this limit 
must involve at least two sufficiently hard photons, arranged such that they recoil, at least 
partially, against each other. Naively, in the classical limit, such configurations are dominated 
by photon emission off both fermions. This motivates why radiation beyond the kink is absent 
in the TV-decay spectra. Along the same lines of reasoning, such double hard photon emissions 
are decreasingly probable with rising lepton masses. However, since in the present state only 
approximated matrix elements up to 0{a) are included in the program these double emissions 
are not described correctly yet. 

In Fig. [T] also different treatments of higher order matrix elements are exhibited: photons 
emitted solely according to the purely soft YFS eikonals (left panel) are contrasted with correc- 
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(a) soft only 



(b) soft & collinear 
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Figure 1: Photon radiation in leptonic decays of Z (upper panel) and W bosons (lower panel) 
for different leptons, including fictional heavy r's in a range of masses. In the left panel, (a), 
C = 1, i.e. photon generation according to the YFS form factor only is depicted, whereas in 
the right panel, (b), corrections up to 0{a) are included using the dipole splitting functions and 
the virtual corrections, cf. Sec. \4.3\ and \4.1\ All distributions are normalised on the total decay 
width of the decay into the respective lepton and lepton-neutrino pair. The infrared cut-off in 
all cases is set to O.lGeV. 



tions due to the approximated matrix elements presented in Sec. 14.11 The former distribution, 
labelled with "soft", thus is correct in the soft limit but it is inadequate for the description of 
hard, collinear photon radiation. This, including virtual corrections of 0{a), is displayed in 
the panel labelled with "soft & collinear". 

The inclusion of these corrections gives reasonably good results as long as most photons are 
soft or if complicated correlations of hard photons are not important. 

5.1.2 Comparison with other codes 

After checking the physical sanity of the implementation in principle, results obtained with 
PHOTONS++ are now to be compared to those from other, established and dedicated Monte 
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Figure 2: The total photon energy in the decay frame in W ^ tv^-. In the left panel, (a), 
the distributions generated by WiNDEC (black), PHOTONS++ (red) and AMEGIC++ (green) are 
depicted, where the latter has been rescaled with the true average photon multiplicity. In the 
right panel, the relative deviations of WiNDEC (black) and PHOTONS++ (red) with respect to the 
rescaled matrix element result of AMEGIC++ are displayed. 



Carlo event generators capable of describing QED effects in the decays of W and Z bosons, 
in particular with the WiNDEC package [Mj. This program aims at the description of the 
production and decay of VT-bosons in hadronic collisions. WiNDEC performs the decay of the 
VT-boson into lepton-neutrino pairs including QED corrections summed in the YFS-approach 
and corrected by exact 0{q) real emission and virtual correction matrix elements. They are 
obtained for the decay only in the narrow width approximation, i.e. only the W ^ iu decay is 
taken into account. Furthermore, in this section, the results of PHOTONS++ are compared with 
the exact, fixed order, one-photon emission results of the SHERPA-inherent matrix element gen- 
erator AMEGIC++. However, comparisons with AMEGIC++ are only sensible when the average 
photon number of the process under consideration is low, i.e. when multiphoton emission gives 
a negligible contribution to the differential cross section. Additionally, it should be stressed 
that AMEGIC++ lacks virtual corrections and therefore comparisons are sensible for normalised 
distributions only. 

The channel best suited for comparing all three generators is W ^ tv^-. Besides the low 
avarage photon multiplicity (with an infrared cut-off of 0.1 GeV multiphoton events make up 
for less than 3% of all radiative events) virtual corrections merely amount to a 1% correction of 
the zeroth order cross section. Furthermore, as discussed earlier, the majority of multiphoton 
events will consist of at most one single hard photon and additional soft ones. Therefore, these 
events will be aproximately described by the hard emission only. 

It should be stressed at this point, however, that there is one fundamental difference in the 
comparison of the various results, related to the way the infrared cut-off is implemented: While 
in WiNDEC the energy cut-off is applied in the rest frame of the decaying W , it is applied in 
the rest-frame of the decaying dipole in both PHOTONS++ and AMEGIC++. 

The distributions generated by all three programs are shown in Fig. O In general terms. 
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the distributions agree reasonably well with each other. There is, however, a slight deviation in 
the region of large radiated energies, where WiNDEC undershoots the results of the two other 
codes on the level of up to 10%. On the other hand, WiNDEC exhibits an overshoot in the very 
low energy bins, for radiated energies around or smaller 5 GeV, which is due to the different 
frames in which the infrared cut-offs are applied. As already mentioned, in WiNDEC this is 
defined in the W rest frame, hence resulting in a fiat hypersurface in this frame. In contrast, 
in PHOTONS++ it is applied in the rest frame of the iy-/-dipole. Subsequently the surface of 
the region cut off by this definition forms a directionally dependent hypersurface in the rest 
frame of the W (observable in Fig. [3] where the cut-off is set to IGeV). The net result is that 
some photons having more than O.lGeV in this frame had less than O.lGeV in the rest frame 
of the dipole, and vice versa. Ultimately, different defintions of the infrared cut-off result in 
different behaviour in the vicinity of this cut-off in nearby frames. The differences are the larger 
the further both frames are apart. On the other hand, the differences at high photon energies 
most likely stem from different mapping procedures in both codes. The mapping procedure 
in PHOTONS++, cf. Sec. 13.3.21 does not involve neutral particles, in this case the neutrino. It 
therefore ensures that the full phase space possible for the radiative decay can be mapped onto 
the leading order one. 

Another feature to study is the dependence of the distributions on the choice of the infrared 
cut-off u. It is employed to separate the divergent region of real soft photon emission, which is 
exponentiated together with the virtual contributions, from the region of the phase space where 
resolvable photons will be generated. Accordingly, this cut-off sets a limit on the number of 
photons to be generated. In Figure[3]the results of this variation on the spectrum of the radiated 
photons' total energy are exhibited for two different final states, electrons (upper panel) and 
r's (lower panel). In the case of the decays W —>■ the two codes, PHOTONS++ (left panel) 
and WiNDEC (right panel) show a similar behaviour: Varying the cutoffs between 1 MeV and 
1 GeV yields stable results in large regions and especially also in the high-energy tail of the 
distribution, whereas differences appear only in the region of small energies, around 1-2 GeV. 
However, in the case of the decays W ^ eu the differences between the two codes are more 
pronounced. Varying the cutoff there yields still comparably stable results for PHOTONS++, 
but the results of WiNDEC show a significant dependence on the cut-off of the order of around 
10%. This is due to the fact that with decreasing fermion mass the effect of the infrared cutoff 
on the avarage photon number increase^. 

In order to choose an optimal value of the infrared cut-off uj there are different considerations 
to be taken into account: On the one hand an efficient generation is desirable, pushing u as high 
as physically sensible, e.g. the detector level energy resolution on soft photons or decay products. 
Along the same lines it should be noted that all photons in the soft (unresolved) region will be 
assumed to yield a negligible combined momentum. Therefore, choosing a comparably large 
infrared cut-off will not have any effect on distributions involving the resolved Bremsstrahlung 
photons, but it will reduce the accuracy of results obtained for e.g. invariant masses of the 
primary decay products. This consideration clearly demands a smaller cut-off. On the other 
hand, when exponentiating the real soft photon emission a factor / '^S{k){e~'y''-l)eitu-k^) 

In fact, this feature was one of the reasons for preferring the r decay channel over the electron channel. 
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Figure 3: Dependence of the total energy of the radiated photons' on the infrared cut-off in 
PHOTONS++ (left panel) and WiNDEC (right panel) for W eue (upper panel) and W — > rUr 
decays (lower panel). The relative difference to uo = {^.1 GeV is shown. 



has been neglected, which is strictly true only for — > 0. Thus, some residual dependence is to 
be expected, even if infrared subtracted matrix element corrections were included to all orders. 
This dependence is of course minimised with small cut-offs. 

5.1.3 Effects of inclusion of exact matrix elements 

Including exact matrix elements, as discussed in Section 14.21 further improves the accuracy 
of the distributions. This is especially true away from the singular limits, where considerable 
differences emerge. This is exemplified in Fig. El where the angular distributions of photons 
in Z ^ decays is depicted. Of course, there is also an effect on the differential decay rate. 
Fig. [6] shows that corrections obtained from the quasi-coUinear approximation, i.e. from the 
approximate matrix element, overestimate the exact matrix element resulting in an increased 
differential decay rate. Even more so in the region of very hard photon emission which, due to 
the angular constraints imposed by the emitter's mass, no longer fulfils the condition {p-k) 0. 
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Figure 4: The total energy of all photons radiated in Z £i. Left panel (a): The same plot 
as in Fig. U^b), but this time the correction is done by using the exact matrix element (solid) 
instead of the approximated one (dashed). The distribution generated using the eikonals only 
is shown as a dotted line. Right panel (b): The relative difference of the distributions obtained 
using the exact and the approximated matrix elements. In both cases again different fermion 
masses have been used. 
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Figure 5: The total energy of all photons radiated inW ^ . Left panel (a): The same plot 
as in Fig. U^d), but this time the correction is done by using the exact matrix element (solid) 
instead of the approximated one (dashed). The distribution generated using the eikonals only 
is shown as a dotted line. Right panel (b): The relative difference of the distributions obtained 
using the exact and the approximated matrix elements. In both cases again different fermion 
masses have been used. 
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Figure 6: Angular distributions of the emitted photons in Z it, using exact and approximated 
matrix elements. In the left panel (a) and the right panel (h), the cases Z ^ ee and Z ^ tt 
are exhibited using the eikonals only (dotted lines) and corrections through exact (solid) and 
approximated matrix elements (dashed). In both plots the leptons sit at 9 = and 9 = n. 



Here the full matrix element exhibits some destructive interference between the two relevant 
diagrams which is of course absent in the treatment through the dipole splitting kernels. This 
leads to a slightly earlier drop-off of the differential decay rate w.r.t. radiated energy than with 
the approximated matrix elements. 

Further, the absence of interference terms in both the eikonals and the quasi-coUinear ap- 
proximation leads to an overestimation of radiation at large angles. Because of only small 
correlations between the energy of the photon radiated and its angular distribution this over- 
estimation leads to an almost constant decline in the differential decay rate w.r.t. the photon 
energy when corrected by the exact matrix element. Of course, while this effect is small in the 
decay channel Z e^e~ , it increases with the mass of the emitter and when a larger fraction 
of the radiation is radiated at large angles. Nevertheless, for very high emitter masses (cf. the 
fictive r with m,- = 40GeV in Fig. Hj) the approximation proves useful again. This is due to the 
dominance of the soft logarithms over the quasi-collinear ones in this limit. 

5.2 Other channels 

Finally, a short overview over other interesting cases is given. In principle, PHOTONS++ can 
handle any possible final state configuration in single particle decays independent of its charge. 
Thus, it is well suited to address all r- and hadron decays, which will be the topic of this 
section. 

5.2.1 J/^ decays to leptons 

First of all, consider the case of J/^ ^ il, which is topologically identical to leptonic Z-decay, 
but nonetheless very important for the calibration of detectors and as a background source of 
leptons. In Fig. [7] the decay channels J/ip ^ e^e~ and J/ip ^ ^~^fi~ are investigated and the 
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Figure 7: The total energy of the radaited photons in the rest frame of the decaying J/ ip vector 
meson for different lepton pairs (electrons in red, muons in green) in the final state. C = 1 
(dotted) is contrasted with C = 1 + PI/Pq, where Pi is calculated in the quasi- collinear approx- 
imation (dashed) and with the complete real emission matrix element (solid). In all cases, the 
distributions are normalised on the width of the inclusive decay into the respective lepton pair, 
and the infrared cut-off has been fixed to uj = IMeV. 



effect of 0{a) corrections is scrutinised. Again, tlie kinematic limit at half the mass of the 
decaying particle produces a visible and prominent kink. Due to the much smaller mass of the 
J/ ip compared to the Z mass, the effects of the higher muon mass are much more pronounced, 
both in the sharpness of the kink and the quality of the quasi-coUinear approximation. 

5.2.2 B D*+ pions and semileptonic B decays 

Another system to demonstrate the versatility of PHOTONS++ are 5-decays because of its 
manyfold topologies in the final state. 

In Figure [H] semi- leptonic decays of B^ mesons into D~ scalars and D*~ vectors are displayed. 
The resulting distributions are similar for e or fi being the lepton. This is because in both cases 
the bulk of the radiation is emitted off the lepton and the amount of phase space open for 
bremsstrahlung is of similar magnitude. Only the avarage photon multiplicity is noticably 
affected by the difference in mass between the electron and the muon. The r-channel on the 
other hand presents itself differently due to the mass of the tau being comparable both to 
the mass of the B^- and the D(*^~-mesons. This not only leads to the near absence of soft 
bremsstrahlung above the infrared cut-off, as compared to the other semi-leptonic channels, it 
also results in a completely different radiation pattern: The bulk of the photons is radiated 
in between both dipole particles and not primarily coUinearly. Furthermore, the absence of 
interference terms in the radiation off the lepton-scalar pair in contrast to the lepton-vector 
pair is plainly visible for both e and /x. Again, the relevance of the interference terms is small 
for the r-mode due to its radiation being dominated by the spin-independent soft terms. The 
exact matrix element correction also shows the shortcomings of the dipole splitting functions 
in this case as they fail to predict the excess of hard radiation for the eletron. This attribute is 
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Figure 8: Semi-leptonic decays D^iu and D*(2010)~£z/ for different leptons and 

with different matrix element corrections. The solid line corresponds to the correction using the 
full matrix elements in the point-like hadron approximation, the dashed line corresponds to the 
dipole splitting kernels neglecting interference terms and the dotted line corresponds to using the 
eikonals only. The angular distributions are shown in the i — D^*^~ rest frame with the lepton 
at 9 = 0. Again, the infrared cut-off was set to IMeV. 



shielded in the muon case by its aheady comparable large mass. However, the total radiative 
decay rate is nearly uneffected by this. The spin-dependence of the dipole approximation is 
also suppressed by the large mass of the D~ and D*~ , respectively, hence the small difference 
of both cases in that approximation. 

As an example for dealing with multiple charged particles in the final state, decays 
into a D* accompanied with various numbers of charged and neutral vr's have been chosen. 
The results are on display in Fig. [HI where the total radiated photon energy and the angular 
distribution of the photons are depicted. The orientation of the final state momenta has been 
chosen in such a way that configurations of the same multipole structure differing only by 
a neutral pion have a similar momentum distribution within the multipole, but still letting 
the tt" have a no n- vanishing effect. The most prominent feature in the distriubtion of the 
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Figure 9: The total photon energy in the rest frame of the decaying meson for different 
numbers of pions in the final state (upper plot) and the angular distribution of this radiation 
in the multipole rest frame with the D*~ at 9 = (lower panel). For the multi-body final 
states the same kinematic configurations have been used, as detailed in the text, to yield easily 
interpretable results. For identical multipoles similar final state momentum configurations with 
non-vanishing vr" momentum have been chosen to increase comparability. The infrared cut-off 
in all cases has been set to uj = IMeV. 



total energy of all photons in the B meson's rest frame is the receding kinematic limit for 
the total energy, it is independent of the momentum layout within the multipole. It is due 
to the decreasing amount of phase space open for bremsstrahlung with an increasing number 
of pions. On the other hand, while the total energy available for the photon decreases, the 
amount of Bremsstrahlung increases with the number of charged particles involved. Switching 
from a dipole (two charged particles) to a quadrupole (four charged particles), the probability 
of double hard photon emission is increased due to favourable momentum configurations among 
the strongly radiating pions. Additionally, since the pions are spin-0, their photon distribution 
is generated exclusively by a product of eikonal factors. Furthermore, the angular distributions 
of the emitted photons are shown. There, the differential cross-section is integrated over energy 
and the azimutal angle. For better interpretability these distributions are plotted in the rest 
frame of the multipole. The D*(2010)~ allways rests at 6* = 0. Due to its large mass, compared 
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Figure 10: The total photon energy in A"*""*" p'^n'^ in the rest frame of the decaying A 
baryon is exhibited. The infrared cut-off was set to IkeV. 



to the pions, very little radiation is emitted in its direction. In contrast, all charged pions are 
plainly visible as peaks in the spectrum. However, their respective mass cones are hidden due to 
the azimutal integration unless the pion sits at 6' = vr, as is the case in the dipole configurations. 

5.2.3 A++ J9+ 7r+ decays 

A rather exotic decay for the purpose of this publication is the decay A"'"^ p+vr"*", due to 
its lack of neutral particles. This case is presented in Fig. [TUl where the total energy and the 
angular distribution of the emitted photons are exhibted. However, this channel leaves only 
very little phase space open for photon radiation. Thus, coUinear enhancement for the p"*" and 
the A"''+ should be negligible. 

5.2.4 T decays 

The leptonic r decays are an example of a final state containing multiple neutral and massless 
particles. This has the effect that the leading order decays do not have a fixed momentum 
distribution among the primary decay products leading to a smearing out of the sharp kink 
at ^mr, as depicted in Figure [TH Because of the relatively small r-mass and the considerable 
fraction of momentum carried by the neutrinos the effects of the different masses of the elec- 
tron and the muon are plainly visible, in the photon energy spectrum as well as the angular 
distribution. 

Furthermore, the branching fraction of radiative leptonic decays in /i and r decays (with at 
least one photon with E.y > lOMeV) has been checked against PDG values [2Sj, cf. Tab. [21 

6 Conclusions and outlook 

In this publication a new implementation of the YFS approach to the description of higher- 
order QED corrections in particle processes has been presented in the form of a Monte Carlo 
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Figure 11: The total photon energy in t~ i~ Vi in the rest frame of the decaying r lepton 
is shown in the left panel. In the right panel the distribution of the photons ' polar angle is shown 
in the t — i rest frame with the t at 9 = 0. In both plots, the solid line shows the distribution 
corrected with the exact matrix element and the dotted line the one using the eikonals only. The 
infrared cut-off was set to IMeV. 







r(T^ei/e!^r7) 








r{T—*fii'fj,i'T,incl.) 


PDG 


0.014(4) 


0.09(1) 


0.021(3) 


PHOTONS++ 


0.0147(1) 


0.0999(3) 


0.0233(2) 



Table 2: A comparison of the branching ratios of the radiative leptonic jj, and t decay mode 
(E^ > lOMeV) in relation to their inclusive leptonic mode calculated by PHOTONS++ and the 
PDG world avarage. The number in brackets reflects the absolute error on the last digit. 



code. It is a part of the multi-purpose event generation framework Sherpa since version 1.1 
and allows for a simulation of photon radiation in particle decays. This is an important effect 
with important experimental consequences. The huge advantage of the YFS approach is that it 
can be systematically improved order-by-order in the electromagnetic coupling constant, such 
that its accuracy can be increased to match exact results at in principle any given perturbative 
order. Thus, in terms of possible accuracy, the YFS approach clearly reaches beyond typical 
parton-shower based algorithms. Some of the effects due the inclusion of exact perturbative 
results have been studied in this publication. 

In contrast to another recent implementation of the YFS approach in SOPHTY, here, in 
PHOTONS++, there is no constraint in the number of particles produced in the decay, i.e. 
PHOTONS++ stretches beyond the level of 1 — >■ 2 decays. This is possible due to a new method 
of reconstructing the kinematics after QED radiation has been added to a core process, thus 
shifting its characteristics a posteriori. The corresponding algorithms have been tested and 
validated in detail through comparison with results from other codes and experimental data. 
Some of the results have also been presented in this paper. 
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It is anticipated that in the progress of the further development of Sherpa also its modules 
will be improved; in the case of PHOTONS++ this will mainly involve the addition of an increasing 
number of exact higher order results. Some of the most relevant 1 — 2, such as generic V — FF 
matrix elements with adjustable couplings, cf. Sec. 14.11 S — >• FF and S — > SS, as well as more 
dedicated W ^ iu, t ^ ii^ei^r, S Siu and S — > Viu are already present. Others will need 
to be added. The structure of the code also permits the inclusion of form-factors to take into 
account the composite nature of hadrons. 
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A The YFS-Form-Factor 

In this appendix, the cancellation of virtual and real soft singularities will explicitly be per- 
formed and the YFS-Form-Factor will be calculated. As already defined in Sections [2] and [3] 
the YFS-Form-Factor Y{Q) reads 

i<j 



where the virtual infrared factor is given by 



B{Pi,Pj) 



57r 



-ZiZjOiOj 



d'^k ( 2p,9i - k 



2pj9j + k 



p \k^ - 2{k ■ p,i)ei k"^ + 2{k ■ pj)ej 



and the real infrared factor reads 



B{p„p„n) = —ZiZ,e,ej d''k6{k'){i-e{k,n)) 



Pi 



{Pi ■ k) {pj ■ k) 



As before, Zi and Zj are the charges of particles i and j in units of the positron charge, 
respectively, and the sign factors 6',^ = ±1 for final (initial) state particles. Again, VL is the 
"unresolved" region of the phase space for the soft photons. In this form the divergences need 
to be regularised, which can be achieved by either introducng a ficititous small photon mass A, 
as in the original YFS paper [8], or through dimensional regularisation. In both cases, however, 
the limited real emission phase space VL will lead to potentially large logarithms. 

After performing the momentum integration, the virtual infrared factor can be written as 



B{Pi,Pj) 



ZiZjOiOj 
2^r 



, rriimj 
In^-^ 



1 /2 1 

+ l{p^■PJm I d^^ + i 

J Px 



dx In 



rriimj 



where 



P'x 



{PiOi - PjOj) + x{pi9i + pjOj 



and 



contains the infrared divergence. Similarly, the real infrared factor reads 

1 -2 1 



B{Pi,Pj: 



UJ 
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, UJ , rriimj 



In ^ 
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dx 
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G(x 

pI 



ln% 

pI 
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with 



Px 



iPi + Pj) + xipi - Pj) 



■ In z — + In 



and cu is the momentum cut-off specifying ft in the frame B is to be evaluated in. Furthermore, 

G{x) 

with 



Wx 



1-Px 



Px^ 



Pi 



{Ei + Ej) + x{Ei - Ej) 



Combining both terms to the YFS-Form-Factor the divergences cancel and a finite result 
is obtained. The remaining parameter integrals do not give rise to further divergences as long 
as > 0, i.e. as long as the emitting particles are massive. Thus, taken together, the YFS 

form factor reads 



Y{pi,pj,u;) = 



a 



TT 



ZiZjOiOj 



In^ 



Wi-Pj) 



+ \{Pr-Pjme,e_ 




dx In 



uiimi 



n^e{x[x'^) 



x'2- x[){pi+pjy 



dx 



\nx^ 

pI 



+ G'(l) + G'(-l)-(p,-p,) /dx^^ 

J Px 



where x'12 are the roots of p'^ with x'-^ < x^- The general case cannot be evaluated in closed 
form. This is due to the fact that the term 

1 



/ 



dx 



pI 



although completely finite, can only be evaluated analytically for the dipole in its rest frame 
or in the rest frame of one of either of its constituent particles. This can only be achieved if 
there is one dipole only. All other cases need to be evaluated numerically. 



A.l Special cases 

A. 1.1 Decay into two particles with [pjOi -\-pj6jY < 

If the multipole consists of only two particles in the final state, e.g. for decays of the type 
Z — > then there is an analytical solution in the rest frame of the dipole formed by the two 
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charged particles. In the high-energy hmit, given by Ei ^ rrii for both QED corrected charged 
particles, the critical term above can be written as 



{pi ■ pj) I da;^^ 



= 6^ • 



Therefore, in this case, the full YFS form factor reads 
Y{p„pj,uj) = 



— ZiZjOiOj 



1-ln^^^^ ln^ + ln|iln^-lln^^ 



rriimj 
^ rriimj 



1 - — 

^ 6 



This result in the high-energy limit agrees with the result stated in [8]. 



A. 1.2 Decay of a charged particle with one charged final state with {piOi+pjOjY = 

A similar, but nonetheless different case occurs for the decay of a charged particle into a final 
state involving only one charged particle, e.g. the case of VT-decays, W iuf. Then, in the 
corresponding dipole's rest frame neither mw Ew nor {pi9i + PjOj)^ < and therefore this 
case is different from the one above. In this case, for {pw — PiY = 0, 



YJu) 



a 

TT 



2 1 - In 



mw 
mi 



This result of course differs from the result in [23] since both results are given in different 
Lorentz-frames. Also, if in this process a photon is radiated, then {pw — PiY = '^{Pu ■ V-t) > 
and the YFS-Form-Factor takes a different a form. 



A. 2 The full YFS form factor 

Here the complete solutions to analytically integrable parameter integrals in the YFS form 
factor are given. In the following, using the invariance of YiVL) under the interchange of 
Pi ^ Pj) the labels pi and pj are chosen such that Ej > Ei. It is useful to define 



p1 -p]^ '^sJiPt-PjY - pip] 



,2^2 
3 

Xl,2 



i.Pt - Pj) 

as the roots of p^ and 



pf -p'j^'^J (Pi -PiY - pIp^ 



{P^+P,f 
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as those of -p'"^ in case of QiQj = +1, satisfying xi^2 ^ [—1)1] ^12 ^ (— 1>1)) respectively. 
It holds that Xi, ^2 > and X2, x[ < if {pi — pj)"^ < and < Xi < X2 and < x[ < x'2 if 
{Pi—PjY > 0. These difference in the relations between xi and X2 necessitate the differentiation 
of distinct cases in the calculations. 

If {pi—PjY — then xi^2 are not defined. If OiOj — —1 then p'^ — p^ and x'12 are meaningless, 
leading to another set of distinct cases. 

When evaluating the first set of the parameter integrals that fact simphfies matters a lot 
resulting in 

1 ..,2 1 



7^e I OiOj 



In '^■ 
dx- 



+ / dx- 



In^ 



A2 



^= 0. 



-1 



Pi 



Otherwise, the evaluation is more complicated and involves shifting the poles at x'^ 2 off the 
real axis. The solution then is 

1 U-2 1 



He I BS, y dx^ + 

STT^e {x\x'^ 
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In 



(4 - A ) {Vi + Pj ) ^ (a^l - X2 ) {pi - Pj ) 



ln|x,|(Li. (^)-Li. (^)) 



- In \x2\ (Li2 ( - 



X2 



- Li2 



X2 + 1 
X2 



In any case, the last piece of the divergence has cancelled, leaving finite terms negligible in the 
high energy limit. 

The other integral containing p'^ is to be evaluated next. In total there are three cases to 
consider. 

• OiOi = +1 



Tie 



da; In 



rriimj 



^ 2\n^ll±^ + \u\{l-x'^){l-x'i)\-x\\n 



l-x[ 


— X2 In 


1-4 


l + x[ 




1 + 4 



4. 



AUthough, there again are poles within the range of integration the integral over them is 
finite. 

• 9i9j = —1. The range of integration does not comprise any poles and, thus, is real, giving 



1 

/ 



dxln 



iriim-i 



= 21n 



+ In [(1 - xl){l - xj)] + Xi In 
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l+Xi 



1 — Xi 



+ X2ln 



I + X2 



1 - X2 



4. 



Evidently, the case {pi —PjY = has to be treated separately. It yields 
1 



/ 



dx In 



mjUii 



2ln— — ^ + ln|l-a;^| +a;„ln 



1 + Xr 



1 Xqn 



-2. 



where Xp = — pi_pi ■ In decay matrix elements it is not kinematically possible to also have 



rrij = nij. 



The last integral that is generally solveable analytically differentiates even more cases. The 
easiest to solve is the case of Ei — Ej, as it is occuring in leptonic Z-decays. Here, E^ is 
independent of x, thus giving 
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{xi - X2){pi-Pjf 
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In 
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For all other dipoles three distinct cases appear: 

• {Pi-PjV < 
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dx- 
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With the definitions for xe and Xp from above it allways holds that xe > Xp > 1, thus 
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The last integral can generally only be solved numerically. This is due to the complexity of 
/3x. If, however, the dipole is in its rest frame or in the rest frame of one of its constituents, 
there are analytical solutions. Because PHOTONS++ allways treats multipoles in their rest 
frames solutions for the integral will only be given in that frame. Two important cases are: 
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with p=\M = M andE = Ei = Ej. 



• Leptonic W-dec&j {rrii <^ rrij = mw) 

1 

/ da:^ ^ A [h-' + Li2 (-2)] • 
J p% 
-1 

B Transforming the phase space elements 

This section details the phase space manipulations neccessary for the implementation of the 
YFS algorithm in form of a computer code. 

B.l Rewriting the phase space element in other frames 

As discussed in Sec. 13. 2[ the phase space integral with the phase space element 

d$ = d^d^ui^T^YKVC+PN-Pc-PN-K) 
n 

= n 

i=l 

has to be transformed to explicitely be in the chosen frame, the multipole rest frame. This can 
be achieved by using the identities 

1 = ^y"d^(pc+Piv)d^Pcdmi,,p5(^^(pc+Piv)) 5((pc+P7v)'-M2) 

X 5' [Pc + Pn- S{{pc + Pc? - ml,^^) 6 ((pc + Vn?) ■ 

and 

1 = [d'x6 - l] 6' (^L-\pc + Pc)) . (78) 

Here, mM,p is the invariant mass of Pm = Pc + Pc and M is the invariant mass of the initial 
state Pc +Pn- As before, pc and pjy and Pc and P/v are the sums of the initial and final state 
charged and neutral particles' momenta. The first identity basically amounts to extending the 
integration to an integration over the full phase space including the initial particles. The second 
identity, taken from [26], involves a Lorentz-transformation, denoted by L~^, being the boost 
into the rest frame of x. Applying this boost on the phase space integral of course is a valid 
operation, since the full expression at this point is formulated in a Lorentz-invariant way. The 
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result of this Lorentz-transformation, after inserting both identities, reads 

d$ = {2n)^d%d<!>k [ d\pc + Pn) d^Pc dml^^ d^x ^\pc + Pn - Pc - Pn - K) 

X 6'' {L{pc+Pn)) S{{pc+Pn)^ - M^) 6^ (Pc + Pn - J^P^ 
X 6{{pc + Pc? - mlj e {{pc + PnT) 

X M ^ - 1 V' (^(Pc + Pc) 

Reordering and using the identity 
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yields 
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The last line can be further simplified by using the identity of Eq. ( ITSll again and by integrating 
over M^. Now, the other integrations can be performed, first over [p + Pn), then over P and 
finally over m\jp. This results in 

d$ = (27r)M$pd$fc^^53(2 5^Pi-P^ + i?-p^)5n5^p, + irj -mM , 

where m?^^ = [pc + Pc)"^ = {2J2Pi ~ Pn + K — pi^Y = Pli the invariant mass of the 
QED-corrected multipole. 
Finally, the identity 

^ + ^)^ - ^t') = 2(Pg + po+JfO) ^'^- - ^'.«) 

will be used, where q = Pq + p% = itlm,p and where all zero- components are taken in the 
rest frame of Pm = Pc + Pc- Therefore, 
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The phase space element d$ has thus been exphcitely rewritten in the rest frame of the multi- 
pole, at the cost of a Jacobian. 

Similarily, the zeroth order uncorrected cross section can be transformed to 

d$o = {27cyd<l>g 6^ (pc+Pn-Qc-Qn) 

where mM,q is the invariant mass of the uncorrected multipole and the and Q^r are taken 
in the Qm rest frame. 



B.2 Rewriting the phase space element in terms of the undressed 
momenta 

In both cases the manipulations can be done in close analogy to the unitary algorithm of 
[27]. The neccessary manipulations are easiest done backwards, starting with the phase space 
integral in terms of the and defining n = + ^at to be the number of final state particles. 



B.2.1 Mixed multipoles 

In this case the starting point reads 
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This can be recast into a better form by inserting the identity 
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and by expressing the (5-function fixing in terms of the kinematically relevant variables 
and p^. This then yields 



5'\Qm)S{Q'm-Q%-p'c) 



= Jdufl d% d'p, 5 {q^ - m^) e(g°) 5' {p, - uq, + R) 5 U - 



x5 



A M2 + lu^qi- ncK j - ^ ^/^j|TCu^^^^ - 
\ \ c J C,N 



\ 



E«: -E' 



c 



C,N 



X 



PcPc 



pm 

Pc c,N Pi . 



Integrating over d^qi and d^^, using 5 {x^ — x^) Q{x) — ^S{x — Xq), and integrating over u 
yields 
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i=i \ ' / ^ 2./4,fp,+ 
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where in the integration over u the second last 5-function of the hne above has been used. 
Furthermore, in this transformation, an identity similar to flHUj) . arising when defining u in 
terms of pj, has been employed. A rearrangement of terms and a suitable transformation of the 
last ^-function in terms of Pm yields 
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Here, the identity 
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has been used. Reversing the procedure allows to express the phase space element through the 
undressed final state momenta as 
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B.2.2 Final state multipoles 

The transformation will be done using the same techniques as above. Starting from 

i=l 

i=l \ C / \ C 

Again, similar identities to ( 1791) and ( IHOi) will be used, but due to the different mapping scheme 
they now read 
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And, as before, the (5-function over Q%j is expressed in the kinematically relevant variables qf. 
This then yields 
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Integrating over d^gj and u yields 
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where, again, the second last 5-function has been used in the integration over u. Additionally, 
an identity similar to (!82|) . arising when defining u in terms of pj, has been used. Rearranging 
terms leads to 
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where the identity 



= V + 



has been used. Reversing the procedure allows to express the phase space element through the 
undressed final state momenta as 

— * 777 

d$ = d$,d$,(27r)^5='(QM)5(Q°M-Qc) 



Piv _ % " r 

Pn'pn _ sr^ pm J-J- 

C Details on the photon generation 

In this section the generation of the photon distribution is detailed. 
C.l Avarage photon multiplicity 

The average photon multiplicity n is the avarage of the Poisson distribution before it is corrected 
by the various weights. It is therefore not immediately connected to the true avarage photon 
multiplicity of the final event. Nonetheless, it is an integral part of the generation procedure. 
An analytical result in closed form is available for both dipoles and multipoles. However, the 
calculations for multipoles are more involved as the integrations do not nicely seperate as they 
do in the dipole case in the chosen frame. Thus, as a starting point the analytical result for 
the dipole in its rest frame will be given. It reads 

n = / -jTrSqik) = ZiZ29i92 m In- 



' TT ^ ^ ^ ^ u^,^ + (l-/3i)(l-/32 

where ujmm is the infrared cut-off and cjmax is the maximal kinematically allowed photon energy. 
The latter can be determined by setting the rescaling parameter u to zero in Eqs. (IHUj) and fHTj) . 
respectively, and by assuming single photon emission. Additionally, = 

In the case of a multipole, the integral over the photon energy can still be separated, as 
long as the soft photon region is sufficiently well-behaved. This is the case, if G(fc, Q) forms an 
isotropic hypersurface in the frame of the integration. However, the angular integration still 
remains to be done: 

—eik,n)s,{k) 



In 5^ V Z^Z,9^9, (Sn- [ dl^-^^-^ 
47r^ t^min \ J [qi ■ ek){qj ■ ek) 
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polar axis 



Figure 12: Sketch of how the axes are chosen in the angular integration in multipoles. 



Choosing different orientations of tlie polar axes for eacli interference term of every con- 
stituent dipole, all angular integrations can be done analytically. Although this may sound like 
quite an ad-hoc procedure, it is completely valid and simplyfies the integration immensely. The 
orientation for each of the interference terms is thus chosen to be such that both momenta lie 
symmetrically in the unit sphere, both forming an angle aij with the polar axis, see Fig. \n\ 
Therefore, by this choice. 



where again is with e\ = 0, cf. Eq. ( l5Ti) . and the further parameters are given by 



{Qi ■ Qj) 

{qi ■ efc) 



EiEj (1 — ajOj + bibj) 

Ei{l — ai sin ip sin 9 — hi cos 9) 

Ej {1 — ttj sin Lf sin 9 + hj cos 9) , 



' sin aij and 6j j = jSij cos aij . 



With these choices the last integral reads 




dip / d^ sin 9 



1 



(1 — Oj sin (y9 sin 9 — hi cos (1 — aj sin ip sin 9 + hj cos 9) 







Using the decomposition 



1 1 



hj{l — Qi sin ipsm9 — hi cos 9) hi (1 — aj sin ipsm9 + hj cos 9) 

{hi — hj) + 2hihj cos 9 



hihj (1 — sin (y9 sin 9 — hi cos 9) {1 — aj sin Lp sin 9 + hj cos 9) 
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and Qibj — ajbi, this can be easily integrated giving 



EiEn 



(Qi ■ ek)(qj ■ efe) 



27r 
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^Cj-Dj + Ej + 



B{2Cj-Dj)-A(Dj-2Ej) 
2^B'^Cj-ABDj+A^Ej 



^B^C, - ABDj + A^E, A-B^c, + D,+E, + ^ 



B{2Cj+Dj)-A{Dj+2Ej) 
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with 



A = 6,-6. 



B 



2bibj 



Ci,j — 1 



D. 



E, 



^2 I 1,2 



Upon closer examination it can be seen that for Oj,- — > the result of the dipole case is recovered. 



C.2 Photon energy 

Due to the dccompostion of the integration over the photon energy and the integration over 
the unit sphere, the photon energy distribution and the photon angular distribution can be 
generated seperately. Of course, this independence of distributions is no longer true after the 
reweighting procedure, but it alleviates the generation of the crude distribution. 

In the imlementation presented here, the photon energy is distributed according to p, 
generated through 

7 / '-^max 

V<^min 

where 7?. is a uniformly distributed random number on the interval [0,1]. 



C.3 Photon angles 

Similar to all other parts of the photon distribution, the photon angles are also generated 
according to Sq{k). For this, the relevant function is recast into the form 

_ y 

i-Pf I 2(1 + A/3,-) 1-/3] 

(l-/3iCose)2 (1- (3i cose) (1 + (3 j cose) (1 + (3j cos e)^ ' 
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where 9 is some polar angle w.r.t. the dipole axis in the dipole rest frame. In this frame, the 
generation of the azimuthal is trivial - it just follows a fiat distribution in [0, 2it]. The polar 
distribution above can be bound from above through the interference term. This allows to 
generate the true distribution by generating the angle according to the interference term and 
applying a hit-or-miss rejection. The interference term can be decomposed analogously to the 
general case above into two independent terms according to 

1 m 



(1- /3iCos^)(l+/3jCos^) + (3j ^icose) A(l + /3jCos^) 

The cosine of the polar angle, cos^^, is then generated to either of the two terms, i.e. it is 
generated according to (1 — A cos 6*)^^ with probability 



Pi = 



In l±§i + In i 
l-3i ^ 1 



and according to (1 + (3j cos 9) ^ with probability Pj = 1 — Pi, selected through a random 
number. These angles can be generated by 



cos^ = 



in the former case and 



cos^ 



Pi 



1 



1 + A 



n 



A- 



in the latter. TZ again is a uniformly distributed random number on [0, 1] . The correction 
weight for obtaining the full distribution reads 



2[l+J,Jj) 



w 



(l-ftcos6>)2 ~ (l-/3iCos6»)(l+/3j COS0) (l+/3j 008 0)2 
2(l+ft/3j) 



< 1. 



(l-/3iCos6»)(l+/33Cos6») 

The azimutal angle (f is distibuted uniformly. 



C.4 Photons from multipoles 

In a multipole configuration again the photons are generated according to Sq{k). The integral 
over photon energies can still be seperated from the angular integrations, decoupling the gen- 
eration of the energy of a single photon as above. However, its angular distribution is very 
complex. But due to 

^qi^) = ^S{qi,qj,k) 

i<j 



50 



the photon angles are distributed according to 

This is nothing but a sum of angular distributions of different dipoles which are not in their 
respective rest frame. 

Subsequently, one of those constituent dipoles is chosen with the probability 

Then, photon angle generation can proceed as above in the rest frame of the dipole. To obtain 
the right distribution in the rest frame of the overall multipole, a null-vector of unit length is 
created in the rest frame of the dipole using the newly generated angles (p ant 9. Then this 
null vector is boosted into the rest frame of the multipole. It now has the angular distribution 
according to its constituent dipole in this frame. Since it is a null vector it has the properties 
of a photon and only needs to be rescaled to the energy generated earlier. 

D Massive dipole splitting functions 

The massive dipole splitting functions are needed for the calculation of the approximation 
to the infrared subtracted single hard photon emission matrix element (31 . They are taken 
directly from [17] for spin-| emitters and are generalised from [19] for all other cases. Problems 
arising during this generalisation are related to the fact that these splitting functions for spin-1 
particles are only given for massless gluons and that all initial states are considered massless 
as well. The extension to radiation off massive spin-1 particles is rather straight forward by 
augmentation with a simple mass term. The extension to massive initial states is less clear 
since decay matrix element are far off the massless initial state limit. However, the decaying 
particle is allways much more massive than its decay products when those are supposed to emit 
hard bremsstrahlung. Thus, photons are predominantly emitted at large angles to the initial 
state resulting in negligible contributions from these splitting functions. Hence, they can safely 
be omitted. 

Also, velocity factors from [19] have been omitted. They were introduced to facilitate the 
analytic integration and change neither the infrared nor the quasi-collinear limit. They only 
result in a different interpolation inbetween. The same is true for the factor Rij in the massive 
fermion splitting function of [T7]. Nonetheless, here this factor is kept because of the direct 
applicability of these splitting functions to the completely massive splitting. 

Three cases need to be differentiated regarding the state, initial or final, the emitter and 
spectator are in. The fourth case where both emitter and spectator are in the initial state lies 
outside the present applicability of this program, it will therefore be omitted. 



rii 



Kj 
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To repeat the notation, pi is the 4-momtentum of the emitter, pj that of the spectator and 
k is the emitted photon. All massive dipole splitting functions will be given, in that order, 
for spin-0, spin-i and spin-1 emitters. Since there are no massive dipole splitting functions 
available for emitters of spin-| or spin-2, their emissions have to be described by the soft limit 
only. Of course, it is allways possible to implement exact process specific matrix elements. 



Final State Emitter, Final State Spectator 



with 



and 



with 



1 



{pi ■ k)Rij{yij] 
1 



2 ^ _ _ mf 

l-Zij{l-yij) {Pi-k)_ 



1- Zij{l- yij) 1- Zkj{l- Vij) 
Pik 



Vij 



PiPj+Pik + Pjk 

Pipj 

PiPj +Pjk 

1 — Zij 



Rij{x) 



Pij = Pi+Pj + k 



(pi + k) ■ Pj 
^(2m| + i^2(l-x))'-4P, 



m] = 2 {piPj + Pik + Pjk) 



wherein the photon is massless, X{x, y, z) is the Kallen-function. 



{Pi ■ k)_ 



Final State Emitter, Initial State Spectator 



gij{Pi,Pj,k) = dij {PhPj^k) 



{Pi ■ k)Xij 

1 

{Pi ' k)Xij 



2 Xij Zij 



^ 

{p^■k). 
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2 Xij Zij 2 Xij Zj^j 



~l~ 2,ZijZi^j 4 



(Pi ■ k). 
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SlS2 


X{pi,si;p;p2,S2]Cl,cr) 


+ + 
+ - 


/iiyU2?7^CL + yU^?7ir/2Ci? + CrS{+;pi,p)S{-;p,P2) 
CLfJ-ivSi+;p,P2) + CRfi2vSi+;pup) 



Table 3: X -Functions for different helicity combinations. Missing combinations can be obtained 
using the simultaneous replacements + ^ — and L ^ R. 



with 



Xij 



PiPj +Pjk-pik 
PiPj +pjk 
PiPj 



PiPj +Pjk 

Initial State Emitter, Final State Spectator 

The emitting particle is allways assumed to be much heavier than its decay products resulting 
in its contributions to the real emission corrections to be negligible. Thus, 

gij{Pi,Pj,k) = gtj''\Pi,Pj:k) 

is set irrespective of the emitter's spin. 



E Basic building blocks For matrix element calculations 

In this Appendix a short summary on the defintions of the basic building blocks (cf. [20121]) for 
the calculations of exact matrix elements will be given. Additionally, techniques to incorporate 
propagators into that scheme will be reviewed. 

X- Function 

The X-function is a contraction over a ferimonic current coupled to a vector with an arbitrary 
structure of the vertex. 

X {pi,si]P]P2,S2]Cl,cr) = u{pi,si) :^[clPl + CrPr]u{p2,S2) , 

where u{pi, Si) may be a particle or anti-particle spinor, = and Pr = The vector 

dotted into the 7-matrix may be a momentum vector or a polarisation vector. For the 
explicite calculation of the X-Function see Table El 
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SlS2 


Y{pi,suP2,S2;cl,cr) 


+ + 
+ - 


CLSi+;puP2) 



Table 4: Y -Functions for different helicity combinations. Missing combinations can be obtained 
using the simultaneous replacements + — and L ^ R. 



5i SoS'iS A 

^ L'^ A'^ 0^ ^ 


Zim , Si ; Do, So: v^i^ Sq; »4, Sa'. ci^, Cd^; c?^, c^) 


+ + ++ 
+ + +- 
+ + -+ 
+ + — 
+ - ++ 
+ - +- 
+ - -+ 
+ 


2 [5(+;p3,Pi)'5(-;p2,P4)c}f4^ + ^ll^J'2V3mcfc^R + ^^3^^AVlV2c]icf] 

'^V2Cr [Si + ■,Pl,P4)^i3cf + S{ + ;PuP3)^l4C^^] 
2^1 C|? [S{-;P3,P2)f^4cf + 5(-;p4,P2)/i34^] 
2 [^( + ;P4,Pl)5(-;P2,P3)4^ci'^ + /il/^2r/3r/4cP4^ + fl3fl4VlV2CR cj^] 

2r/44^ [5(+;pi,p3)/i242 + 5(+;p2,P3)yUicP] 


-2 [/ii/i4?72?73cPcf + IJ,2fJ'3VmCRC^j^ - /ii/i3?72r/4cPc|^ - IJ,2lJ'4VlV3CRcf] 
"^VsCr [S{ + ■,P4,P2)^^lCf + S{ + ;PuP4)H2Cr] 



Table 5: Z -Functions for different helicity combinations. Missing combinations can be obtained 
using the simultaneous replacements + ^ — and L R. 

y-Function 

The y-function is the pendant of the X-function when the fermionic current is coupling to a 
scalar rather than a vector. 

y {Pi,si;p2,S2]Cl,cr) = u{pi,si)[clPl + CrPr]u{p2,S2) . 
Its explicit calculation is shown in Table HI 

Function 

The Z-function is a contraction over two ferionic currents connected by a massless gauge boson 
(cf. Table El). 

Z {pi, si; p2, S2; P3, S3; p4, S4; c]^ , c^j^; cf, cf) 

= U{pi, Si)7^ [c^Pl + 42 Pr] U{p2, S2MP3, 53)7^ [cfPi + C^Pr] u{p4, S4) . 

5- Function 

For the calculation of the above spinoral products it is useful to define the S'-Function 

Sis;pi,p2) = u{pi,s)u{p2,-s) . 
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Its two possible forms for given pi and p2 are 

, X _ ^ {Pi ■ ko){p2 ■ ki) - {pi ■ h){p2 ■ ko) - ieai3y5PiP2kokf 
'-'{-riPl:P2) — ^ 

e/ N _ o ■ ^o)iP2 ■ ki) - {pi ■ ki){p2 ■ ko) + ieaP'y5PiP2k2kf 

^[~iPl^P2 — , 

where ko is an arbitrary null vector (/cq = 0) and ki satisfies the relations kf = —1 and 
{ko ■ ki) = 0. Furthermore, 



Vi 

It is also useful to define the quantity 



y/2{pi ■ ko) . 



rrii 

IJ'i = ± — 

Vi 



where ± refers to particles/anti-particles. 
Fermionic Propagators 

These propagators can be incorporated using the following identity: 

1 ± —= u{p, s)u{p, s)+ [It —f= I v{p, s)v{p, s) 



= IE 

s 

This allows to cut the line and replace it with a sum of external particles 
Bosonic Propagators 

Bosonic propagators can be incorporated by writing out their Lorentz-structure explicitely. This 
is trivial in Feynman gauge, if the vector is massless . Massive propagators are best included 
in unitary gauge, since then no additional goldstone boson exchange has to be included. 
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